i/i 


1D-A192  476 8  NUMERICAL  TREATMENT  OE  THE  FLUID/ELASTIC  INTEEACE 

UNDER  RANGE-DEPENDENT  <U)  VALE  UNIV  NEU  HAVEN  CT  DEPT 
OF  COMPUTER  SCIENCE  E  C  SHANG  ET  AL  JAN  88 
UNCLASSIFIED  VALEU/DCS/RR-SSS  F/G  8/il  NL 


rn.BCB.COP)  f  £ 


^  ET  VERil^ 


ELECTE 

MAR  2  5  B88 


A  Numerical  Treatment  of  the  Fluid/Elastic  Interface 
under  Range  —  dependent  Environments 

E.  C.  $hang1and  Ding  Lee2 
Research  Report  YALEU/DCS/RR-558 
January  1988 


DISTHSUTION  STATEMENT  A 


Approved  fox  public  talsoset 
Distribution  Unlimited 


YALE  UNIVERSITY 

DEPARTMENT  OF  COMPUTER  SCIENCE 


88  8  18 


ABSTRACT 


f?:u 

\  N 


i 


i 


& 


The  technique  introduced  by  McDaniel-Lee  for  the  handling  of  the  fluid/fluid  inter¬ 
face  boundary  under  range-dependent  environments  is  extended  to  handle  the  horizontal 
fluid/elastic  interface  boundary.  Representative  wave  equations  of  the  parabolic  type  are 
considered  in  both  fluid  and  elastic  media.  The  required  interface  conditions,  (1)  continuity 
of  vertical  components  of  displacement,  (2)  continuity  of  vertical  components  of  stress,  and 
(3)  horizontal  components  of  stress  vanish  on  the  interface,  are  satisfied  with  this  numerical 
treatment.  A  complete  theoretical  development  is  presented  along  with  a  test  example  to 
demonstrate  its  validity. 


*  A  action  F  o* 


I  NT l£  CRAM 


'  OTiC  TA9 

I  U: 

!  >.  .  „  ...... 


.1!  0. 


I  Cy  fiOy  l4f _ 

A\ ,.ii id 

(”  "  ;  Avc  1  if..! ;  or 

!  v  .  c.j| 


A  Numerical  Treatment  of  the  Fluid/Elastic  Interface 


under  Range  —  dependent  Environments 


E.  C.  Shang'and  Ding  Lee3 
Research  Report  YALEU/DCS/RR-558 


January  1988 


DTIC 

/'TjNELL-CTU'r 
VX  f/Ah  2  5  19881* 


1  !'  f  V"-l 


‘Present  Address:  NOAA/ERL  Wave  Propagation  Lab.,  Boulder,  CO  80303. 
’Present  Address:  Naval  Underwater  Systems  Center,  New  London, CT  06320. 


I.  INTRODUCTION 


Modern  computational  techniques  Improved  much  of  the  efficiency  of 
range-dependent  wave  propagation  models.  The  efficiency  as  well  as 
usefulness  of  the  range-dependent  model  can  be  further  enhanced  If  a 
capability  can  be  incorporated  to  handle  the  fluid/elastic  interface.  This 
paper  introduces  a  numerical  treatment  of  the  horizontal  fluid/elastic 
Interface  boundary.  Since  the  numerical  treatment  introduced  by 
McDaniel-Lee  [1]  to  handle  the  fluid/fluid  Interface  in  1982,  Interest  was 
on  the  rise  in  searching- for  efficient  methods  to  handle  the  fluid/elastic 
Interface.  Some  contributions  were  made  In  the  seismology  science  to  treat 
the  fluid/elastic  Interface.  R.  Stephen  [2]  developed  finite  difference 
methods  to  treat  this  problem  dealing  with  a  hyperbolic  wave  equation. 

J.  T.  Kuo  and  Y.  C.  Teng  [3]  applied  the  finite  element  as  well  as  finite 
difference  schemes  extensively  to  solve  the  same  problem  as  above  but 
dealing  with  an  elliptic  wave  equation.  The  above  techniques,  though 
workable,  are  not  simple  to  adapt  Into  any  existing  range-dependent  model 
without  requiring  excessive  efforts;  moreover,  these  techniques  are  by  no 
means  simple.  The  numerical  treatment  we  Introduced  In  this  paper  Is  based 
on  the  standard  Parabolic  Equation  (PE)  in  the  fluid  medium  Introduced  by 
Tappert  [4]  and  on  the  coupled  parabolic  equations  In  the  elastic  medium 
derived  by  McCoy  [5].  The  fluid/elastic  Interface  requires  three  conditions 
to  be  satisfied  on  the  Interface,  l.e.,  (1)  the  continuity  of  vertical 
components  of  displacement,  (2)  the  continuity  of  vertical  components  of 


stress,  and  (3)  horizontal  components  of  stress  must  vanish  on  the 
Interface.  These  conditions  were  derived  to  be  consistent  with  the  PE 
representation  In  both  fluid  and  elastic  media.  McDanlel-Lee's  technique  was 
modified  to  treat  these  conditions  numerically.  This  modification  allows  the 
existing  impllclte  finite  difference  ( 1 FD )  [6]  marching  scheme  to  be  applied 
systematically.  A  test  problem  with  known  solution,  given  by  Ewing  and 
Press  [7],  Is  used  to  examine  the  validity  of  this  development. 


II.  BACKGROUND  SUMMARY 

Since  the  new  treatment  to  the  fluid/elastic  interface  boundary  is  an 
extension  of  the  McDaniel-Lee  technique,  It  is  desirable  that  the 
McDaniel-lee's  treatment  to  the  fluid/fluid  Interface  boundary  be  briefly 
reviewed. 

We  use  u(r,z)  to  Indicate  the  wave  field,  the  pressure,  in  a  2- 
dlmenslonal  medium,  depth  and  range.  Thus,  u(r,z)  satisfies  the  parabolic 
wave  equation 

ur  *  a(ko.r,z)u  +  b(ko,r,z)  uzz  ,  (2 

where  a(k0,r,z)  *  1kQ(n2(r,z)  -  l)/2 

and 


b(k  ,r,z)  -  1 /( 2k  ) , 
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where  kQ  Is  the  reference  wavenumber,  n(r,z)  stands  for  the  index  of 

refraction  which  is  defined  as  a  ratio  of  the  reference  sound  speed  c  and 

o 

the  sound  speed  c(r,z). 

At  the  ocean  bottom,  the  change  of  sound  speed  and  density  form  an 
Interface  (see  Figure  1).  From  one  medium 
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Figure  I:  Interface  Boundary 

to  the  next,  at  each  Interface,  the  interface  conditions  must  be  satisfied, 
i.e.,  the  pressure  and  normal  components  of  particle  velocity  are  continuous 
at  the  Interface. 

The  standard  PE,  Eq.  (2.1),  does  not  contain  the  density.  In  order  to 
satisfy  the  interface  conditions,  McDaniel-lee  developed  a  special  equation 
with  density  variations  to  represent  the  wave  field  on  the  Interface.  It 
turned  out  that  this  special  equation  is  again  a  PE. 


In  developing  this  Interface  PE,  McOaniel-Lee  applied  the  Taylor 
series  expansion  to  the  points  near  the  Interface  and  then  match  the  fields 


m 


on  the  Interface  which  Is  denoted  by  z0.  A  clear  configuration  Is  given 
In  Figure  2. 
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Figure  2:  The  Interface  Between  Two  Media 

In  carrying  out  the  matching  process,  let  us  describe  the  Interface 
mathematically  below. 

The  continuity  of  pressure  requires  that 


u(r,z0)  -  u2(r.zB), 


the  continuity  of  normal  component  of  particle  velocity  requires  that 


(2.2) 
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Note  that  the  Interface  Is  assumed  horizontal, 


In  every  medium  the  field  u(r,z)  satisfies  the  PE,  (2.1).  Therefore, 

In  medium  1,  u^r.z)  satisfies 

(»Vr  *=  a1  (k0,r,z)  u1  +  b1(k0,r,z)  (u,)^  ,  (2.4) 

where  a1(kQ,r,z)  «  1kQ(n2(r,z)  -  1 )/2  •  (2.5) 

and 


b1(ko,r,z)  =  1/(2kQ).  (2.6) 

Using  the  first  three  terms  of  the  Taylor  expansion  for  (u^)^ 

upon  (u  )  and  solve  for  (u. )  ,  then  substitute  them  into  Eq. 
i  m  i  zz 

(2.4),  we  find 

3lJi  l  u  i 

az-  "  2b^  (Vr  '  2b^  alul  *  h  (U1  ^Vm-l*  (2‘7) 

where  h  =  Az. 

Similarly,  In  medium  2  use  a  three-term  Taylor  expansion  for 

^u2^rrH-l  upon  ^u2^m  an<*  ^°^ow  *he  same  procedures  as  carried  out  in 
medium  1 ,  we  find 


_h_ 
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<U2>r+ 


JL 

2b, 


a2U2  +  h  ((u2)m+l 


u2) . 


(2.8) 
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The  first  interface  condition  (2.2)  allows  one  to  write  ■  u2  -  u  on 
the  Interface.  Then,  multiply  both  sides  of  (2.7)  by  p2  and  multiply  both 
sides  of  (2.8)  by  p1 ;  then  the  second  Interface  condition  (2.3)  allows  the 
above  results  to  be  equal.  After  simplification,  the  McDaniel-lee 
horizontal  Interface  wave  field  is  obtained  between  two  fluid  media  to  be 


Note  that  the  density  Is  assumed  to  remain  constant  in  each  medium. 

III.  FIELD  REPRESENTATIONS  IN  FLUID/ELASTIC  MEDIUM 

Dealing  with  the  fluid/elastic  Interface,  two  media  are  Involved, 
i.e.,  a  fluid  medium  and  an  elastic  medium.  The  field  representation  In  the 
fluid  medium  by  the  parabolic  approximation  is  a  scalar  PE  while  in  the 
elastic  medium  the  representative  wave  equation  Is  a  vector  PE.  A  standard 
derivation  of  the  scalar  PE  in  fluid  medium  can  be  found  in  references  4 
[Tappert]  and  8  [ Lee-Si egmann ] .  The  vector  PE  has  been  derived  by  a  few 
authors,  though  their  derivations,  and  method  of  derivation  differ  from  one 
another.  These  derivations  are 


(1)  Lander  and  Claerbout's  equation  [9],  derived  using  dilatation 


and  rotation, 

(2)  Hudson's  equation  [10],  derived  using  displacement,  and 

(3)  McCoy's  equation  [5],  derived  using  dilatlonal  and  shear 
potentials. 

In  this  paper  we  use  potential  representation  to  deal  with  the 
fluid/elastic  Interface,  thus,  it  Is  appropriate  that  we  adopt  the  vector  PE 
developed  by  McCoy. 

In  the  fluid  medium  the  parabolic  wave  equation  was  expressed  in  the 
2-dlmensional  cylindrical  coordinates  where  we  use  ^(r.z)  to  Indicate  the 
potential;  in  the  elastic  medium  we  use  $2(r,z)  and  4'2(r,z)  to  represent 
the  potentials  there.  In  the  above  notation,  the  subscript  "1"  Is  used  to 
indicate  the  fluid  medium,  and  the  subscript  "2'*  to  Indicate  the  elastic 
medium.  We  use  to  Indicate  the  displacement  In  the  fluid  medium  and 

0?  to  Indicate  the  displacement  In  the  elastic  medium.  Their 
relationships  with  the  potential  are  given  by 


D]  (u1,w1)  =  grad  ^(r.z), 

(3.1) 

°2  (u2,w2)  =  9rad  *2  (r,Z)  +  rot  *Vr,Z)' 

(3.2) 

Equivalently  we  can  write  the  above  as 
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(3.4) 


r* 

El 

iV“ 

$ 

1 

$ 


I 

I 

I 


>: 


3*2  3*2 

ar  "  az 


3t)  2  3  ^2  l 

dz  ar  r  ^2 


In  the  fluid  medium,  the  potential  ^(r.z)  satisfies  the  Helmholtz 
equation  below. 


V2*](r,z)  +  k2  (r.z)^  =  0. 


(3-5) 


Following  the  derivation  of  the  PE  in  the  fluid  medium  by  Lee-Siegmann, 
the  wave  field  *.j(r,z)  obeys  the  following  decomposition 


(1), 


«(r.z)  =  A^r.z)  H^'(k 


‘or)  -  Vr'2>  4). 


(3.6)a 


where  kQ  Is  the  reference  wavenumber  and  (kQr)  is  the  zeroth 

order  Hankel  function  of  the  first  kind.  In  the  foregoing  development,  all 

Hankel  functions,  in  their  asymptotic  expansion  form,  have  a  common 

[7  '  1  7 

multiplicative  constant  *1-  e  ;  this  can  be  Ignored  for  simple 

calculation.  This  treatment  allows  Eq.  (3.6)a  to  be  written  as 

j— —  ikQr 

*(r,z)  =  A(r,zj^|~k  r  e  (3.6)b 
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Assuming  that  A^r.z)  Is  weakly  dependent  in  range,  the 
representative  PE  in  the  fluid  medium  is  obtained,  by  the  parabolic 
approximation 


V,  *  *,  ^ 


12  2 
where  a.,  =  ££  (kl  (r,z)  "  ko) 
0 


(3.7) 


(3.8) 


(3.9) 


In  an  inhomogeneous  elastic  medium,  4>2(r»z)  and  no  longer 

satisfy  the  two  independent  Helmholtz  equations  because  of  the  coupling 
effects.  Following  the  derivation  of  McCoy,  the  propagation  of 
time-harmonic  stress  waves  through  an  Inhomogeneous,  linear  elastic  solid 
medium  which  is  locally  Isotropic,  Is  governed  by  the  equation 


-  -  2  -  2 
(X2  +  y  •  0?  +  Uy  v  +  pw  D?  +  F?  =  0 


2  M2V  2 


'2  2 


(3.10) 


where  \2,  v2  ar®  Lame  parameters  and  p 2  is  the  density  such  that 


X2  "  X2  +  cr(r,z)i’ 


v2  *  v2  ^  +  cw^r,z^* 
p2  -  [1  +  cp(r,z)], 


(3.11) 

(3.12) 

(3.13) 


WiTf; 


where  e  Is  a  unit  vector. 


In  the  above  equation,  the  upper  bar  indicates  the  spatial  average  and 
c  with  a  subscript  indicates  the  perturbation  with  respect  to  that  subscript. 


Define 


-  2 
p2  03 


X2  +  2^2 


and 


assuming  that 


(3.15) 


(3.16) 


4>2(r,z) 


and 


+2(r,z) 


A2(r,z)  H</>(kD  r)  ~  A2(r,z) 


B2(r,z)  (k$r)  ~  B2(r,z) 


(3.17) 


(3.18) 
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IV.  FLUID/ELASTIC  INTERFACE  CONDITIONS 


In  this  section,  we  derive  the  fluid/elastic  interface  conditions 
associated  with  the  representative  PE's.  Furthermore,  from  these  interface 
conditions,  a  system  of  three  equations  which  relate  the  fluid  potential 
A1  to  the  elastic  potentials  A?  and  B2  and  their  derivatives  on  the 
fluid/elastic  interface  will  be  derived. 


We  begin  by  discussing  the  fluid/elastic  interface  conditions  using 
the  expression  of  displacement.  These  conditions,  in  cylindrical 
coordinates,  are 


1.  Continuity  of  vertical  components  of  displacement: 


W1  =  W2‘ 


2.  Continuity  of  vertical  components  of  stress: 


(4.2) 
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3.  The  horizontal  components  of  stress  must  vanish  on  the  interface: 


hi  + 

\az  ar  J 


(4.3) 


In  terms  of  potentials,  the  equivalent  interface  conditions  to  (4.1), 


(4.2),  and  (4.3)  are 


1.  Continuity  of  vertical  components  of  displacement: 


a«l*1  a<t>2  a*2  1 


if  +  FT  +  ?V 


(4.4) 


2.  Continuity  of  vertical  components  of  stress: 


(3\  ,  1  !i_  ,  l3\  1 1  !*>  iS\  ,  5 

V2  r  ar  az2  /'  2  \  ar2  r  4r  az2/  “2  \  az2  4r4z  r  42  / 


(4.5) 


3.  The  horizontal  component  of  stress  must  vanish  on  the  interface: 


Lfh  ih  fVi  ft 

"2  V  4r4z  '  az2  ar2  F  4r 


(4.6) 


Making  use  of  Eq.  (3.6)a  for*^  (3.17)  for  *2>  (3.18)  for  *2,  and  the 
PE's  (3.7),  (3.19),  and  (3.20),  the  corresponding  Interface  conditions  for 


the  PE's  are  obtained  below. 


1.  Continuity  of  vertical  components  of  displacement: 


3A  IT"  3A,  lAnr  ST”  /2B-  \  iA  r 

ir-Vf  if'  “ 


(4.7) 


2.  Continuity  of  vertical  components  of  stress: 


-p-jM  A^  = 


(3  A-  3A,  ,  \  3  A,  k  iAnr 


(4.8) 


(3  8,  3B,\  k  1 A  r 
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3.  The  horizontal  components  of  stress  must  vanish  on  the  interface: 


(4.9) 


where  An*  k_  -  k  ,  A  *  k  -k. 
D  D  0  S  S  0 


From  this  point  on,  since  we  deal  with  the  potentials  and  their  partial 
derivatives  on  the  Interface  boundary  and  at  the  present  range  level,  for 
economy  In  writing  we  drop  the  superscript  n  and  the  subscript  j,  e.g.,  A? 
means  (A^)"  unless  otherwise  specified. 


Following  McDaniel-lee' s  technique  in  fluid  medium,  we  use  the  Taylor 

32A, 

expansion  for  (A.).  .  upon  (A.).,  solve  for  - r  and  substituting  into 

1  j-i  1  J 


■‘Cstvflyfr  Vvfr/.v'fr. ^  ^  '■ 


Eq.  (3.7),  we  find 
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aA1  h_  ^1  h_  *  1  ,A  ,4  v  X 

az  *  2b]  ar  '  2b.,  a1Al  h  (A1 


(4.13) 


where  h,  is  defined  as  the  depth  increment  az.  Later  we  will  use  k  to 
represent  the  range  increment  ar. 


If  we  substitute  (4.13)  into  (4.7),  we  obtain 


L!!i.!!_aA  +  1  (A 
2b]  ar  2b]  a1Al  h  'A1 


dA„ 


"  (A.),  ,)  -  K 


rj-i;  d  az 


+  ks(tF  *  'E,e2)' 


(4.14) 


where 


1  anr 
o  „  D 


(4.15) 


and 


Ik  ia  r 
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(4.16) 


Using  the  finite  difference  for  the  partial  derivatives  in  (4.14),  we  obtain 


2b]  k  ^1  '  Al^“  2b,  alAl  +  h  (A1  ^Vj-l*  “  *0  h  ^ 

+  Ks(k  ^  -  B2]  +  1  k's  b2)- 


)  "  K0  h  l(Vj+l  A2 


(4.17) 
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A  simplification  of  (4.17)  gives 


Pll  Ai  +  Pi,  A"  +  Pi,  8 


12  2 


'13  “2 


RHS  1 


(4.18) 


where 


h  1 


11  2b]  k 


(4.19) 


Pi 2  *  0, 


(4.20) 


Pi  3  =  ~Ks/lc» 


(4.21) 


^2b1  (k  +  V  '  hj  Ai  +  h(Al > j-1  +  h  (( 


<Vj+l 


+Ks  <iks  -V  B2* 


(4.22) 


The  next  two  equations  to  be  derived  Into  the  system  are  based  on  Eqs. 
(4.8)  and  (4.9).  In  those  two  equations,  two  terms  are  Involved,  namely. 


a2A2  a2B2 

iraz  and  ariz  *  We  f1rst  trV  t0  develop  explicit  expressions  for  these  two 


partial  derivatives. 


Applying  the  finite  difference  to  the  partial  derivatives  In  Eq.  (3.20) 


gives 


(B2>jtl  -  <VJ+1  +ka2  (B2>j+l  +kbife)j+l 


(1  +  k  a2)  (B2)J+1  +  k  b 


2  V  azyj+l  +  k  3Z/j+1 


(4.23) 
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»raz  ar  laz  /  ar  h  ((B2)j+l  B2)  h  V  k  k 


-k{ 
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-  ((B2>r  -  b2)| 
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(4.24) 


-  ((Vi*1  *  B*)j 

Applying  the  finite  difference  to  the  partial  derivatives  In  Eq. 
(3.19),  we  obtain 


(A2)j+l  “  (Vj+l  +  k  a2(A2}j+l+  k  b2^  azl)j+l  +  k  cz(  az)j+l  (4-25) 
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Using  the  finite  difference  for  all  partial  derivatives  In  (4.28),  we 


obtain 
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A  simplification  of  (4.29)  gives 
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■p1w  1  '  h2  f  2  kVj+2 


Using  the  finite  difference  for  all  partial  derivatives  in  (4.35),  we 


obtain 


k  Wj+1  k  b2  k2  |(A2)j+2  "  2(A2Jj+l  +  A2 


hk  |k  a2(A2}j+l  k  b2  ^2 

((A2>  -  A^j+  21kD  ((A2>j+1  "  A2) 


+  k  *  I<B2>^1 


"  ^  UB2)j+2  -  2(B2)j+1  ♦  B2]  -  21ks  1  [B2+1  -  B2] 


(ik$)  B2)  Kq£. 


A  simplification  of  (4.37)  gives 


.n+1  .n+1  _n+l 

P31  A1  +  p32  A2  +  P33  B2  “  RHS3 


where 


P31  ■  Q. 
p32  =  —2/ ( hk ) , 

P33  -  2iks  (1/k)  K0S, 


and 


RHS3  "  '  £  b2  <Vj+2 


+  2lKn  h  (A,) 


D  h  '  2' j+1 


jhk  |ka2  “  2k  b2 

'  ('hk  k  b2  “2  hk  +  2ik0  h)A2  +  72  K0S(B 


2;j+2 


(4. 


(4. 


(A. 

(4. 

(4. 


(4.42) 


-  {k  ^  +  Kos|  <B 


2Jj+1 


+<  ^  +  21k. 


k  ‘  iks  ^  K0S  +  k  h  |  B2 ‘ 


Equations  (4.18),  (4.30),  and  (4.38)  constitute  a  system  of  3  equations 


among  which  unknowns  a"+1  ,  A^2,  and  Bp+1  are  to  be  solved. 


From  this  system  we  see  that  a"+  only  appears  in  Eq.  (4.18). 

Actually,  we  need  only  to  solve  a  system  of  2  equations,  i.e.,  Eqs.  (4.30) 


and  (4.38).  This  system  can  be  written  as 


P22  P23  a"+1 


P32  p33  J  |_B2+1 


(4.43) 


After  system  (4.43)  is  solved,  we  substitute  B,  into  Eq.  (4.18) 


to  obtain  a"  by  means  of 


(RHS  1  -  p13  )/  pn 


(4.44) 


To  examine  the  existence  and  uniqueness  of  the  solution  of  system 


(4.43),  we  evaluate  the  determinent  of  (4.43), 


determinant  of  (4.43)  = 


[x2  2iko 


k)  (2iks  k  K0S^ 


,¥? 


,  2  2 

4  -L-  Ks  w  0 

Irk^  ks 


because  KQS  is  complex  while  other  quantities  are  real;  therefore,  the 
determinent  is  not  singular  and  the  solution  exists  and  is  unique. 


V.  A  NEW  COMPUTATIONAL  APPROACH 

Following  McDaniel-Lee  [1],  a  uniform  partition  in  the  z-direction  is 
assumed.  The  h  =  az  is  used  for  the  depth  increment  and  the  integer  index  j 
is  the  interface  boundary.  As  before,  the  superscript  indicates  the  range 
level,  and  the  subscript  indicates  the  depth  level.  It  is  also  understood 
that  if  both  the  superscript  and  the  subscript  are  dropped,  it  denotes  the 
field  at  (nAr,  jaz),  i.e.,  A  «  Aj. 


Our  approach  can  be  described  by  the  following  statement: 

Solving  the  representative  parabolic  wave  equations  in  the 
different  media  by  means  of  an  implicit  finite  difference  (IFD) 
marching  scheme,  applying  the  field  values  on  the  interface  as 
boundary  information  from  fluid  and  elastic  media. 


t 


We  proceed  to  discuss  the  meaning  of  the  above  statement.  Solving  the 


representative  parabolic  wave  equation  using  an  implicit  finite  difference 


scheme  in  a  marching  process  requires  the  initial  field  values  at  range  r 


plus  the  surface  and  bottom  boundary  information.  The  surface  boundary 
point  at  the  present  level  is  denoted  by  (A^  and  at  the  advanced 
level  by  (A-j)^  the  bottom  boundary  point  at  the  present  level  is 
denoted  by  ( A1 ) j  and  at  the  advanced  level  by  (A-j)j*1.  The  IFD  scheme 
predicts  the  wave  field  at  the  advanced  level,  r  +  a r,  regardless  whether 


the  medium  is  fluid  or  elastic.  If  the  medium  is  fluid,  there  is  only  one 
PE  to  solve;  if  the  medium  Is  elastic,  there  Is  a  system  of  2  PE's  to  solve 
in  this  case  the  field  Is  a  vector  containing  components  A2  and  B^,  each 
component  is  a  subvector.  Dealing  with  the  solution  of  the  entire  problem, 
the  surface  remains  unchanged,  but  an  Interface  boundary  comes  into 
existence.  This  Interface  boundary  separates  a  fluid  medium  and  an  elastic 
medium.  One  crosses  the  interface  boundary,  the  density  and  the  sound 
speeds  change.  In  the  elastic  medium  two  sound  speeds  occur,  the  speed  of 
P-wave  c_  and  the  speed  of  S-wave  cc.  Initial  field  values  at  r  are 
used  along  with  surface  points  ( A-j ) ^ ,  ( A^ ) and  Interface  boundary  points 
( A-j ) j ,  ( A-j ) j to  predict  the  wave  field  (rQ  +  Ar)  in  the  fluid  medium. 

The  (A2)j.  (A2)j+1,  are  used  as  surface  points  along  with 

initial  field  values  at  the  same  range  level.  Beyond  interface  boundary, 

the  boundary  points  ^2)[ottm.  <A?Citom.  d2>bottM»'  "^bottom 
are  generated  artificially.  This  setup  allows  the  same  IFD  procedure  to 

solve  a  system  of  PE's  in  the  elastic  medium  in  the  same  manner  as  in  the 


fluid  medium.  The  key  treatment  Is  to  determine  the  interface  boundary 
values  ( A^  > j •  ( A^ ) j ,  and  (B2)j  which  are  related  by  the  system 


(4.43)  and  Eq.  (4.44).  With  these  interface  boundary  values,  the  IFD  can 

inarch  In  range  to  predict  the  wave  field  in  the  fluid  at  the  next  range. 

This  Is  where  the  McDaniel-Lee  Interface  treatment  Is  extended  to  handle  the 

fluid/elastic  horizontal  Interface.  Note  that  (A.).  is  used  as  a  bottom 

'  J 

boundary  point  to  solve  the  PE  In  the  fluid  medium  while  ( j  an<* 

(B2)j  are  used  as  two  "surface"  points  to  solve  the  system  of  two 
parabolic  equations  in  the  elastic  medium. 


For  better  understanding,  the  diagram  below  explains  our  description  and 
will  help  clear  up  our  early  statements. 


m 


Let  (A*,)".  (A*,)j+1  be  2  vectors  associated  with  range  levels  n 

and  n+1  respectively.  These  2  vectors  contain  only  the  nonzero  components 

at  the  interface  boundaries,  i.e.,  ( A1 ) ^  and  (A1)j+1. 

Then  the  matrix  representation  for  the  fluid/elastic  interface  problem 
can  be  expressed  by 


where  (^j  an<*  (^j  are  ^  vectors  having  the  same  structure  as 
( )  j  except  the  nonzero  elements  are  the  first  components  and  D,  E,  F 
are  tri-diagonal  matrices,  and  E1 ,  are  sparse  matrices  whose  non-zero 
elements  appear  on  the  diagonal  and  lower  diagonal. 

Note  that  the  points  which  influence  the  fluid/elastic  interface 
boundary  are  (A.^,  (A^,  and  ^B2^j  at  a11  ran9e  levels.  It  is 
very  clear  that  their  relationships  are  defined  by  the  system  (4.43)  and  the 
Eq.  (4.44).  Once  these  boundary  values  are  obtained,  the  IFD  scheme  can 
march  forward  in  range. 
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VI.  A  NUMERICAL  EXAMPLE 


This  section  examines  the  validity  of  our  approach  and  an  example  whose 
solutions  [6]  are  known  is  used  as  a  demonstration. 

First  of  all,  in  the  case  where  the  elastic  waves  are  absent  which 
implies  that  is  zero.  Under  such  an  environment  only  one  PE  is  needed 
to  represent  the  wave  propagation  In  the  fluid.  Thus,  the  equation  (2.1)  is 
the  representative  equation.  Furthermore,  Eq.  (4.6)  becomes  zero 
Identically  on  both  sides.  Eq.  (4.4)  reduces  to 


3A]  3A2 

az  "  az 

and  Eq.  (4.5)  reduces  to 


(6.1! 


p2A1 


p/2 


(6.2! 


Eqs.  (6.1)  and  (6.2)  are  equivalent  interface  conditions  to  (2.2)  and  (2.3) 
in  the  fluid  medium.  Then  all  Taylor  expansions  can  be  applied  following 
the  McDaniel-Lee  technique. 


Next  we  use  an  example  discussed  by  Ewing  et  al  [7]  to  computationally 
demonstrate  the  validity.  This  example  neglects  the  branch  line  Integrals 


and  Is  feasible  for  large  values  of  range;  thus,  It  is  very  suitable  for  our 
PE's  (far-fleld).  The  solutions  $lt  4*2  are  9^ven  f°r  Helmholtz 
equations.  For  consistency  we  derive  the  corresponding  solutions  to  the  PE. 


The  exact  mathematical  expressions  for$.j,  $2,  and  4»2  are 


4> 


1 


1(«t-k  r-J) 

e  *-j (kn)  s1n(  end)  sin(  £nz) 


0  <  Z  <  H  = 


d  =  source  depth 


(6.3) 


*2 


1(wt”*cnr“4) 


*2(kn)s1n(  end)e-n(z_H) 


z  >  H 


(6.4) 


1  (ut 


V'!> 


-  (,( Z-H) 

*?(fcn)sin(  i  d)e 


*"d  $n«*n  i0--1  • 

vi 


(6.10) 


nn  "  kn  V  ~  ~2 


(6.11) 


S,  '  kn  V 1  "  c2/B22  * 


(6.12) 


where  n  is  a  subscript  which  indicates  that  the  quantity  is  to  be  evaluated 

at  k  *  k  ,  where  k  are  the  roots  of 
n  n 


p,  * 


1  “-r  ?  tanUH)  -  4k2  ;-(2k2  -  ) 

o  n  H  5  A  £ 


For  computational  simplicity,  we  use  one  mode,  i.e.,  n  «  1.  After 
separating  the  H^1 ^  (k*r)  (k*  *  kQ,  ks>  or  kp)  and  the 
time-harmonic,  the  corresponding  PE's  A1 ,  A2»  and  B2  to  Eqs.  (6.3), 
(6.4),  and  (6.5)  become 


(6.13) 


Ik  1(k  -k  )r 

S  *  "  0  Vk,> 

.  2.45”  p1(kn"  k°)r  «  (k  1  sinf  A  dle'n(z_H) 

*2  =  H~^  e  •2(kl)  s1n(  tnd)e 

Ik  1  ( k  -k  )r  .  , _  ... 

B2  =  2wV^7  6  "  S  W  Sln{  *nd)e”  ^  H 


*2(k1)  sin(  6nd)e 


Hk- -k  )r 


B2  -  2w\jf  e 


♦2(ki >  s1n(  *nd>e"  ;(Z"H) 


(6.14) 


(6.15! 


(6.16 
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has  the  same  definition  as  (6.9)  but  has  the  subscript  n  ■  1 


Next  select  the  congressional  wave  velocity  v-j  in  a  fluid  of  1500  m/s 

2 

with  compressibility  *  v^  and  a  water  density  of  *  1.0.  c  is 
the  phase  velocity  such  that 


f  -  c  lcn/2.. 


We  follow  the  case  I  of  Press-Ewlng's  to  select  o^,  02*  and  vi  such 


<*2  >  02  >  c  >  v] 


We  make  the  following  choices: 


f  =  68.03  Hz 


H  *  Zj  *  100m 


d  =  z$  *  25m 


v1  -  1500m/s 


Wv 


(A^  Is  calculated  by  means  of  formula  (4.44)  where  RHS1  is  defined  by 
formula  (4.22),  p13  is  defined  by  formula  (4.21),  is  defined  by  formula 

(4.19),  and  B^+1  Is  calculated  by  formula  (6.16). 

As  Illustrated  in  the  previous  section,  our  main  effort  is  to  determine 
(A-j)0*1  and  use  it  as  a  fluid/elastic  boundary  Information.  (A.j)j+1 
Is  obtained  In  such  a  way  that  it  is  related  to  the  information  of  elastic 
potentials  (A^)j+1  and  on  the  fluid/elastic  interface. 

System  (4.43)  was  developed  to  relate  these  points.  In  this  test  example,  we  use 
an  accurate  ( ^  and  ( ) ^  at  every  range  from  a  known  solution  as  the 
accurate  boundary  Interface  values.  These  accurate  values  are  applied  at  every 
step  when  solving  system  (4.43).  Results  are  tabulated  describing  the  comparison 
of  computed  field  values  against  the  known  solution  in  dB.  Accuracy  was  carried 
up  to  2  significant  digits  using  the  VAX  11/780  computer. 


(-0.366E-02.-0. 590E-02) 

( -0. 366E-02.-0. 590E-02) 

(-0.669E-02.-0.198E-02) 

(-0.662E-02.-0.209E-02) 

( 0 . 480E— 02 , —0 . 785E— 02 ) 

( 0 . 482E  — 02 , —0 .77  7  E— 02 ) 

( -0. 874E -02,-0. 270E-02) 
(-0.872E-02.-0.276E-02) 

( —0 . 276E— 02 , — C . 428E— 02 ) 
(-0.269E-02.-0.433E-02) 

(-0.464E-02.-0.1 52E-02) 

( -0.487E-02 ,-0.1 54E-02) 

Computed 

Exact 


From  this  selected  test  example.  It  Is  clearly  seen  that  the  numerical 
results,  produced  by  this  model,  agree  satisfactorily  with  the  exact 
solution.  The  results  not  only  demonstrate  the  validity  of  this  model,  but 
also  show  the  correct  computational  procedure  following  the  IFO  procedure. 
This  also  serves  as  an  early  Indication  that  this  model  can  be  readily 
Incorporated  into  the  IFO  code. 


VII.  CONCLUSIONS 


A  mathematical  model  has  been  developed  by  means  of  the  parabolic 


approximation  method  for  handling  the  fluid/elastic  interface.  The  complete 
mathematical  development  plus  the  numerical  example  proved  the  validity  of 


u 

is 


the  model.  As  It  stands  now,  even  though  this  model  is  accurate,  it  is 
limited  to  narrow  angle  propagation  only.  However,  an  Important  feature  of 
this  model  is  that  it  can  handle  a  range-dependent  index  of  refraction  in 
the  elastic  medium.  Moreover,  another  attractive  feature  is  that  this  model 
Is  readily  adaptable  into  the  existing  IFD  code  without  requiring  excessive 
effort. 
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ADDENDUM 


SUBROUTINE  UFIELD 

************************************************************* 

***  USER  STARTING  FIELD 

***  USER  WRITES  THIS  SUBROUTINE  IF  GAUSSIAN  FIELD  NOT  DESIRED 
***  UFIELD  IS  CALLED  IF  INPUT  PARAMETER  ISF  IS  NOT  ZERO 

***  UFIELD  SUBROUTINE  SUPPLIES: 

U  -  COMPLEX  STARTING  FIELD 

★A*********************************************************** 


PARAMETER  MXLYR=101 ,MXN  = 1 0000 , MXSVP= 1 0 1 ,MXTRK=101 ,NIU=1 , 

C  NOU=2 , NPU=6 

COMPLEX  ACOFX , ACOFY , BCOF , BOTX , BOTY , BTA , HNK , HNKL , SURX , SURY , TEMP , 
C  U ,  X ,  Y 

REALM  K0,  K1  ,  KS,  KD,  MU2 ,  LAMNDA2 

COMPLEX  PHII(MXN),  PHI2(MXN),  PSI2(MXN),  CARG,  B1  ,  A1 , 

T1  2 ,  T2 1 ,  T22 ,  T32,  T33,  CHI,  A1JM, 

KDD,  KSS ,  CARH,  CARI 

/USTFLD/  K0,  K1 ,  KS ,  KD,  SIGMA1 ,  CPHI1,  CPHI 2 ,  CPS 1 2 , 

B1 ,  DELTAS,  DELTAD,  A1JM 
EQUIVALENCE  (H,  ZLYR( 1 ) ) ,  (D,  ZS),  (Z,  ZI ) 


COMMON  / I FDCOM/ACOFX, ACOFY, ALPHA, BCOF, BETA (MXLYR) , BOTX, BOTY, 
BTA(MXN) , CO ,CSVP(MXSVP) ,DR,DR1 , DZ , FRQ , I HNK , ISF, ITYPEB, 

I TYPES , I XS VP (MXLYR) , KSVP ,  N,  N1 , NLYR , NS VP , NWSVP , R 1 2 ( MXN ) , RA , 
RHO ( MXLYR ) ,RSVP, SURX, SURY, THETA, TRACK (MXTRK, 2) ,U(MXN) , 

X (MXN ) , XKO , Y(MXN ) , ZA , ZLYR( MXLYR ), ZP, ZS , ZSVP(MXS VP) 

DATA  PI /3. 141 592654/, DEG/57. 29578/ 

C 

C  PUT  THESE  VALUES  IN  TEMPORARILY 
C 

DATA  C/1507.5/,  VI/1500.0/,  ALPHA2/1 725 . 00/ ,  BETA2/1 530 . 00/ , 

D  RHO1/1.0/,  RH02/1 . 97/ ,  IND/0/ 

C 

C 

SEC(ANG)  =  1.0  /  COS(ANG) 

C 

C  ***  STARTING  FIELD  -  EWING  &  PRESS 

C 

AR  =  RA 

IF  (IND  .NE.  0)  AR  =  RA  -  DR 

OMEGA  =  2.0  *  PI  *  FRQ 

K0  =  OMEGA  /  CO 

K1  =  OMEGA  /  C 

KS  =  OMEGA  /  BETA2 

KD  =  OMEGA  /  ALPHA2 

C  WRITE  ( NPU ,  1)  ’K0:  ' ,K0,'K1:  ' ,K1,'KS:  ',KS,'KD:  ' , KD 
1  FORMAT  (2X,A4,E12.6,3X,A4,E12.6,3X,A4,E12.6,3X,A4,E12.6) 


CARG 

=  CMPLX( 0 . 0 , 

( K 1 

-  K0) 

*  AR) 

CARH 

=  CMPLX( 0.0, 

( K1 

-  KD) 

*  AR) 

CARI 

=  CMPLX( 0.0, 

( K 1 

-  KS) 

*  AR) 

B1  = 

CMPLX (0.0,  1 

.0  / 

(2.0 

*  K0)  ) 

C  WRITE  (NPU,  2)  'CARG:  ’ ,CARG,’B1:  ',B1 

2  FORMAT  (2X,A6, ' ( ' , El  2 . 6 , 2X , El  2 . 6 , ' ) ' ,5X,A4, ' ( ' , El  2 . 6 , 2X, El  2 . 6 , 

F  '  )  '  ) 

DELTAS  *  (KS  -  K0)  *  AR 
DELTAD  3  (KD  -  K0 )  *  AR 
DELTADS  =  (KD  -  KS)  *  AR 

C  WRITE  (NPU,  3)  'DELTAS:  ', DELTAS DELTAD :  ', DELTAD, ' DELTADS :  ' 


cc  ccc  cc  cc  non  ta  m  co  co  m 


W  DELTADS 

FORMAT  (2X,A8,E12.6,3X,A8,E12.6,3X,A9,E12.6) 

MU 2  =  RH02  /  BETA2  **  2 

LAMNDA2  =  RH02  *  ( ALPHA2  **  2  -  2.0  *  BETA2  **  2) 
WRITE  ( NPU ,  4)  ' MU2 :  ' ,MU2 , ' LAMNDA2 :  ' , LAMNDA2 
FORMAT  (2X,A5,E12.6,3X,A9,E12.6) 


ARG1  =  C  **  2  /  VI  **  2  -  1.0 
ARG2  =  1  . 0  -  C  *  *  2  /  ALPHA2  **  2 
ARG3  =  1.0  -  C  **  2  /  BETA2  **  2 
ARG4  =2.0-C**2/  BETA2  **  2 

WRITE  (NPU, 5)  ' ARG1 :  ' , ARG1 , ' ARG2 :  ' , ARG2 , * ARG3 :  ' , ARG3 , 

W  ' ARG4 :  ' , ARG4 

FORMAT  (2X,A6,E12.6,3X,A6,E12.6,3X,A6,E12.6,3X,A6,E12.6) 

SIGMA1  =  K1  *  SQRT  ( ARG 1 ) 

ETA1  =  K1  *  SQRT  ( ARG2 ) 

ZETA1  =  K1  *  SQRT  ( ARG 3 ) 

WRITE  (NPU, 6)  'SIGMA1:  ' , S IGMA1 , ' ETA  1 :  ' , ETA1 , ' ZETA1 :  ’,ZETA1 
FORMAT  ( 2X , A8 , E 1 2 . 6 , 3X , A6 , E 1 2 . 6 , 3X , A7 , E 1 2 . 6 ) 

BRACE  =  (RHOI  /  RH02 )  *  (C  **  4  /  BETA2  **  4) 

*  (SIN  (SIGMA1  *  H)  /  (SQRT  ( ARG1 )  *  SQRT  ( ARG2 ) ) 

*  (1.0  +  ARG 2  /  ARG1 )  -  ( <  K1  *  H  *  SQRT  ( ARG2 ) )  /  ARG1 

*  SEC  (SIGMA1  *  H)))  -  4.0  *  (SQRT  ( ARG 3 )  /  SQRT  ( ARG2 ) 

+  SQRT  (ARG2)  /  SQRT  ( ARG 3 )  +  2.0  *  SQRT  ( ARG2 ) 

*  SQRT  (ARG3)  -  2.0  *  ARG4 )  *  COS  (SIGMA1  *  H) 

ZPHI1  =  -((RHOI  /  RH02 )  *  (C  **4  /  BETA2  **  4)  *  ( ETA1  /  SIGMA1 ) 

*  K1  *  H)  /  (1.0  *  BRACE  *  COS( SIGMA1  *  H)) 

CPHI2  =  -((RHOI  /  RH02 )  *  (C  **  2  /  BETA2  **  2)  *  ARG 4  *  K1  *  H) 

/  ( SQRT (ARG 1 )  *  BRACE) 

CPS  1 2  =  -((RHOI  /  RH02 )  *  (C  **  2  /  BETA2  **  2)  *  ( ETA1  /  SIGMA1 ) ) 
/  BRACE 

WRITE  (NPU, 7)  ’BRACE:  ' , BRACE ,' CPHI 1 :  *,CPHI1 
FORMAT  (2X,A7,E12.6,5X,A7,E12.6) 


DO  10  I = 1 , N 
ZI  =  I*DZ 


IF  ((Z  .GE.  0.0)  .OR.  (Z  .LE.  H) )  THEN 

PHII(I)  =  ((2.0  *  PI)  /  H)  *  CEXP(CARG)  *  CPHI 1 

*  S I N ( S I GMA 1  *  D)  *  S I N { S I GMA 1  *  Z) 

*  CSQRT ( CMPLX ( K0  /  K1 ,  0.0)) 


END  IF 

IF  (Z  .GE.  H)  THEN 

PH  1 2 ( I )  =  ((2.0  *  PI)  /  H)  *  CEXP(CARH)  *  CPHI 2 

*  S I N ( S I GMA 1  *  D)  *  EXP ( -ETA1  *  (Z  -  H)) 

*  CSQRT ( CMPLX (KD  /  K1 ,  0.0)) 

PSI 2 ( I )  =  (2.0  *  PI)  *  CEXP ( CAR I )  *  CPSI 2 

*CMPLX( 0 . 0 , -K1 ) 

*  S I N ( S I GMA 1  *  D)  *  EXP ( -ZETA1  *  ( Z  -  H ) ) 

*  CSQRT ( CMPLX (KS  /  K1,  0.0)) 


END  IF 

U ( I )  =  PHII(I) 

CONTINUE 

PHI 2 ( N+ 1 )  -  (  (2.0*PI)/H)*CEXP(CARH)*CPHI2 

*SIN(SIGMA1 *D) *EXP( -ETA1 * ( (N+1 )*DZ-H) ) 
♦CSQRT ( CMPLX ( KD/K 1,0.0)) 


PS  I  2 ( N+ 1 )  =  ( 2 . 0  *P I ) *CEXP (CARI ) *CPSI 2  *CMPLX (  0 . 0 , -K1 ) 

J  *SIN(SIGMA1*D)*EXP(-ZETA1*( (N+1 )*DZ-H) ) 

J  *CSQRT(CMPLX(KS/K1 ,0.0) ) 

KDD  =  CSQRT ( CMPLX ( K0  /  KD,  0.0))  *  CEXP ( CMPLX ( 0 . 0 ,  DELTAD) ) 
KSS  =  CSQRT ( CMPLX (KO  /  KS ,  0.0))  *  CEXP ( CMPLX ( 0 . 0 (  DELTAS)) 

T3 3  =  CMPLX (0.0,  DR  /  (KO  *  DZ  **  2)) 

WRITE  ( NPU ,  8)  ' T33 :  ' ,T33 

FORMAT  ( 2X , A5 , ' ( ' ,E12.6,2X,E12.6, ' ) ' ) 

T32  =  ( CMPLX (1.0,  0.0)  -  T33)  *  PHII(N) 

T32  =  T32  +  T33  *  PHI1 (N-1 ) 

WRITE  (NPU,  8)  'T32:  ' ,T32 
T12  =  CMPLX (0.0,  KS) 

T12  =  KSS  *  T33  *  DZ  *  PSI2(N)  *  (T12  -  CMPLXd.O  /  DR,  0.0)) 
WRITE  (NPU,  8)  ' T1 2 :  * ,T12 

CHI  =  T12  -  T32  /  T33 
WRITE  (NPU,  8)  'CHI:  ' , CHI 

T21  =  KSS  *  T33  *  CEXP(CMPLX( 0 . 0 ,  (K1  -  KS )  *  DR)) 

*  PS  1 2 ( N )  *  DZ  /  DR 

IF  (IND  .EQ.  0)  T21  =  KSS  *  T33  *  PSI2(N)  *  DZ  /  DR 

WRITE  (NPU,  8)  ' T2 1 :  ',T21 

T22  =  KDD  *  T33  *  (PHI2(N+1)  -  PHI2(N)) 

WRITE  (NPU,  8)  'T22:  ' ,T22 

FORMAT  (2X,A6, ' ( ' ,E12.6,2X,E12.6,  '  )  '  ) 

A1  =  T32  +  T22  +  T21  +  T12 

IF  ( IND  .GT.  0)  GO  TO  15 

U  ( N )  =  A1 

IND  =  1 

RETURN 

A1JM  =  A1 

RETURN 

END 


nnnn 


SUBROUTINE  BOON 

***************************************************************** 

***  USER  PREPARED  BOTTOM  CONDITION  SUBROUTINE 

BCON  IS  CALLED  IF  INPUT  PARAMETER  ITYPEB  =  1 
SEE  MAIN  PROGRAM  FOR  DEFINITIONS 

***************************************************************** 

***  SUBROUTINE  RETURNS: 

BOTY , BOTX 

***************************************************************** 


PARAMETER  MXLYR= 1 0 1 , MXN= 1 0000 , MXSVP* 1 0 1 ,MXTRK=101 fNIU=1 , 

C  NOU=2 , NPU=6 

COMPLEX  ACOFX , ACOFY , BCOF , BOTX , BOTY , BTA , HNK , HNKL , SURX , SURY , TEMP , 

C  U ,  X ,  Y 

REALM  K0,  K1  ,  KS ,  KD 

COMPLEX  CARG,  B1 ,  CKS,  CKD,  PI  1  ,  P13,  A1JN,  A1JM,  A2JN,  A2JP1N, 

C  B2JN,  B2JNP1 ,  RHS1 ,  CARH ,  CARI 

COMMON  /USTFLD/  K0 ,  K1 ,  KS ,  KD,  SIGMA1 ,  CPHI1 ,  CPHI 2 ,  CPS  I  2 ,  B1 , 

C  DELTAS,  DELTAD,  A1JM 

COMMON  / I FDCOM/ ACOFX , ACOFY , ALPHA , BCOF , BETA ( MXLYR ) , BOTX , BOTY , 
BTA(MXN) ,C0 ,CSVP(MXSVP) , DR, DR1 , DZ , FRQ, I HNK , ISF, ITYPEB, 

I TYPES , IXSVP(MXLYR) ,KSVP,N,N1 , NLYR , NSVP , NWSVP , R1 2 ( MXN ) , RA , 
RHO ( MXLYR ) , RS VP , SURX , SURY , THETA , TRACK ( MXTRK , 2 ) , U ( MXN ) , 
X(MXN) , XKO , Y ( MXN ) , ZA, ZLYR( MXLYR) , ZP , ZS , ZSVP (MXSVP ) 
EQUIVALENCE  (H,  ZLYR(I)),  (D,  ZS),  ( Z J ,  ZLYR( 1 ) ) 

DATA  P 1/3. 141 592654/, DEG/57.29578/ 

I F ( THETA )  50,100,150 

***  theta  LESS  THAN  0.0.  BOTTOM  SLOPES  UP. 

CONTINUE 
BOTY=U ( N ) 

BOTX= . 

RETURN 

***  TheTA  *  0.0.  BOTTOM  IS  FLAT. 

CONTINUE 
BOTY  «  U ( N ) 

CALL  UF I  ELD 
BOTX  =  A1JM 
RETURN 

***  THETA  GREATER  THAN  0.0,  BOTTOM  SLOPES  DOWN. 

CONTINUE 

BOTY= . 

BOTX= . 

RETURN 

END 


